% multilayer model of temporal spectrum
function [W_aj,params] = temporal_multilayer(params,Z)
% params: parameters
% V: matrix comes from SVD
% Z: matrix for saving j, n and m

lambda = params.lambda;
W_V = params.W_V;
Cn2 = params.Cn2;
L = params.L;
D = params.D;
t_miu = params.t_miu;
f_y = params.f_y;
layerNum = params.layerNum;
L0 = params.L0; % Outer scale
l0 = params.l0; % Inner scale
sp = params.sp;

% Cn2 and wind speed for the multilayer spectra analysis
% dimensions: DeltaZq [m], Cn2q[m^(1/3)], Vq[m/s],zq[m]
MLM = zeros(layerNum,4);

H = (0:layerNum)*L/layerNum;
H_bot = H(1:end-1);
H_top = H(2:end);
MLM(:,1) = (H_top-H_bot)'; % 分块区间的区间大小 [m]
for i = 1:layerNum
    H_mid = H_bot(i) + (H_top(i)-H_bot(i))/2;
    MLM(i,4) = H_mid; %z_q
    if (strcmp(class(Cn2),'double'))
        MLM(i,2) = Cn2;
        MLM(i,3) = W_V;
    else
        intFun1 = @(h) Cn2(h*10^(-3));
        intFun2 = @(h) W_V(h*10^(-3));
        MLM(i,2) = 1/MLM(i,1)*integral(intFun1,H_bot(i),H_top(i)); %Cn2_q
        MLM(i,3) = 1/MLM(i,1)*integral(intFun2,H_bot(i),H_top(i)); %V_q
    end
end
params.MLM = MLM;

% multilayer model for r0m and r0 parameter
r0m = zeros(1,layerNum);
for i = 1:layerNum
    r0m(i) = r0_generate(Cn2,lambda,H_bot(i),H_top(i),sp,L);
end
params.r0m = r0m;
r0 = r0m_merge(r0m);
params.r0 = r0;

% computing the temporal spectrum for each layer of trubulence path
W_aj = 0;
intFun3 = @(h) (1-h/L).^(sp*8/3);
for j = 1:layerNum
    [t_miu2, f_y2] = meshgrid(t_miu/MLM(j,3)*(1-MLM(j,4)/L)^sp, f_y);
    QM = FourierZnkPlns(Z, t_miu2, f_y2, D);
    fun_cal = (abs(QM).^(2)).*(((t_miu2).^2+(f_y2).^2+L0^(-2)).^(-11/6))...
            .*exp(-(2*pi*l0/5.92)^2*((t_miu2).^2+(f_y2).^2));
    cal_integral = sum(fun_cal,1)*(1/D);
    F_aj = 0.00969*(2*pi/lambda)^2*cal_integral;
    W_aj = (MLM(j,1)*MLM(j,2)/MLM(j,3))*F_aj*intFun3(MLM(j,4)) + W_aj;
    fprintf('multilayer calc %3g / %3g \n', j, layerNum);
end

% end of the function
end